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AN  ANALYTICAL  INTEGRATION  OF  THE 
AVERAGED  EQUATIONS  OF  VARIATION  DUE  TO 
SUN-MOON  PERTURBATIONS  AND  ITS 
APPLICATION 

C.  C.  Chao 

The  perturbed  variations  of  the  motion  of  earth 
satellites  due  to  the  sun  and  the  moon  are  derived 
from  a  singly  averaged  disturbing  function.  A 
first-order  solution  is  obtained  by  analytically 
integrating  the  equations  of  variation  including  J^, 

J^,  J3,  and  J^.  The  literal  expansions  are  car¬ 
ried  out  by  computer  in  terms  of  classical  ele¬ 
ments.  The  secular  part  of  the  first-order  solu¬ 
tion  is  included  in  the  reference  orbit.  The  orbits 
of  the  sun  and  the  moon  are  assumed  circular, 
and  the  motion  of  the  moon  is  converted  to  the 
earth  equatorial  system  with  certain  approxima¬ 
tions. 

Results  based  on  the  GPS  {Global  Positioning  Sys¬ 
tem)  satellites  compare  favorably  with  numerical 
integration  for  time  spans  of  up  to  three  years. 

An  algorithm  applying  the  first-order  solution  has 
been  developed  to  achieve  the  desired  strategy  of 
orbit  maintenance  for  the  GPS  Phase  III  system. 

The  analytical  solutions  provide  insight  into  the 
long-term  (10-yr)  variations  of  the  orbit  elements 
of  GPS  satellites. 

INTRODUCTION 

The  need  to  rapidly  predict  the  long-period  and  secular  perturbations  of 
artificial  satellite  orbits  increases  as  more  high-altitude  satellites  with 
long  life  spans,  such  as  the  GPS  (Global  Positioning  System)  satellites, 
are  launched  into  orbit.  By  the  late  1980s,  a  total  of  18  GPS  satellites 
will  be  orbiting  around  the  earth  for  12-hr  periods.  The  ground  tracks 
covered  by  these  satellites  are  required  to  repeat  every  24  hr  to  insure 
a  continuous  global  coverage.  Regular  orbit  maneuvers  must  be  per¬ 
formed  to  maintain  the  specified  repeating  coverage  against  long-term 
perturbations  due  to  the  sun-moon  attractions  and  the  earth  gravitational 
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potential.  The  desiied  strategy  for  orbit  maintenance  maneuvers,  as 
requested  by  the  GPS  program  office,  is  to  minimize  the  total  number 
of  such  maneuvers  during  the  life  span  of  these  satellites.  To  achieve 
this  strategy,  an  analytical  solution  which  predicts  the  long-term  varia¬ 
tions  of  the  orbit  parameters  is  highly  desirable  in  determining  the 
proper  4v  at  the  time  of  each  maneuver. 


During  the  past  two  decades,  efforts  have  been  made  by  many  research¬ 
ers  to  develop  general  perturbation  theories  for  the  long-term  varia¬ 
tions  of  planetary  and  earth  orbiters.  To  name  a  few,  Kozai*,  Kaula^, 
Musen5,  Murphy"*,  and  Estes^  have  either  developed  or  applied  analy¬ 
tical  solutions  to  the  problem  of  third-body  perturbation  for  artificial 
satellites.  Kaufman^,  Uphoff^,  and  Cefola,  et  al.  ®,  have  applied  the 
method  of  averaging,  in  a  numerical  approach,  to  the  prediction  of  the 
long-term  variations  of  planetary  and  lunar  orbiters.  Their  results 
show  a  tremendous  improvement  in  speed  over  that  obtained  from  the 
n-bodv  numerical  integration,  yet  they  do  not  give  the  analytical  ex¬ 
pressions  which  are  preferred  in  this  application. 


The  purpose  of  this  analysis  is  to  further  expand  the  singly  averaged 
disturbing  function  of  third-body  perturbation  arrived  at  by  Kaufman  in 
terms  of  classical  elements,  and  to  apply  the  concept  of  the  interme¬ 
diate  reference  orbit  for  obtaining  a  first-order  solution  by  analytical 
integration.  The  averaged  equations  of  variation  for  zonal  harmonics 
up  to  the  fourth  order  are  included.  The  literal  expansions,  partial 
differentiations,  and  integrations  are  carried  out  by  a  computerized 
series  expansion  technique  employed  by  the  author  in  studying  the  mo¬ 
tion  of  the  Galilean  satellites^.  The  package  of  series  expansion  soft¬ 
ware  used  in  this  analysis  is  the  latest  version  designed  for  the  CDC 
7600  computers  by  Broucke1®.  Although  similar  work  has  been  done  by 
Estes,  the  author  considers  it  necessary  to  carry  out  the  expansion  and 
integration  of  this  particular  application  for  two  reasons.  With  the  help 
of  computerized  expansion  it  would  take  less  effort  and  be  more  con¬ 
venient  to  start  from  the  averaged  disturbing  function  developed  by 
Kaufman  than  to  use  the  analytical  solutions  published  by  Estes.  Fur¬ 
thermore,  Estes  uses  Brown's  lunar  theory  for  the  motion  of  the  moon 
which  carries  a  great  number  of  terms  and  is  deemed  to  be  unnecessary 
for  this  application.  In  this  study,  an  earth  equatorial  coordinate 
system  is  used  with  orbit  elements  of  the  moon  projected  into  this  sys¬ 
tem  using  spherical  trigonometry.  For  time  spans  of  less  than  three 
years,  a  mean  inclination  and  nodal  rate  interpolated  from  the  true 
variation  have  been  found  to  be  good  approximations  for  the  analytical 
integration. 

The  rest  of  this  paper  is  divided  into  two  parts.  Part  I  describes  the 
theoretical  formulation  of  the  first-order  solution  and  compares  the 
long-term  variations  of  the  orbit  elements  of  the  GPS  satellites  pre¬ 
dicted  by  the  analytical  solution  with  numerical  integration.  Part  II 
presents  an  application  of  the  first-order  theory  for  the  desired  stra¬ 
tegy  of  the  orbit  maintenance  of  the  GPS  satellites.  Results  of  the 
investigation  of  long-term  orbit  perturbations  of  the  GPS  Phase  III  sys¬ 
tem  are  discussed. 
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PART  I:  THEORETICAL  FORMULATION 


Averaged  Variational  Equations 


The  disturbing  function  due  to  the  third-body  perturbation  may  be  given 
in  the  earth-centered  coordinates  as 


(1) 


with 

2 

H  -  gravitational  constant  (k  m') 

r*  =  distance  to  the  third  body  (sun/moon) 

r  =  distance  to  the  satellite 

S  =  angle  between  the  two  position  vectors,  r  and  r' 

For  the  sun-moon  perturbations  of  artificial  earth  satellites,  we  note 
that  r/r'  «  1.  Therefore,  we  may  expand  the  square  root  term  in  the 
above  equation  and  neglect  all  terms  of  order  (r/r')3  and  higher. 
Kaufman^  has  carried  out  the  expansion  and  applied  the  method  of 
averaging  to  eliminate  the  short-period  terms.  The  resulting  disturb¬ 
ing  function  is  considerably  simpler. 


with 

A  =  P  •  u ' 
B  =  Q  .  u' 


where  a,  e,  P,  and  Q  are  orbit  parameters  of  the  perturbed  body  and 
a',  n',  r',  and  u1  are  parameters  associated  with  the  perturbing  (third) 
body.  The  u'  is  the  unit  vector  of  the  position  of  the  third  body. 

Kaufman  numerically  integrated  the  variational  equations  of  classical 
elements  derived  from  the  above  disturbing  function  in  terms  of  A  and 
B;  he  then  demonstrated  that  his  method  had  a  500:1  increase  of  speed 
over  an  n-body  numerical  integration  scheme.  Although  his  method  is 
fast  and  reasonably  accurate,  it  does  not  give  the  analytical  expres¬ 
sions  which  are  important  for  analyzing  and  identifying  the  long -period 
terms  of  the  perturbations. 

The  purpose  of  this  paper  is  to  further  expand  the  disturbing  function  in 
terms  of  classical  elements;  thus,  the  variational  equations  (or  per¬ 
turbed  equations)  can  be  obtained  by  partial  differentiation.  With  the 
assistance  of  a  computerized  series  expansion  technique,  we  can  expand 
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the  parameters  (A^  +  B^)  ami  (A^  -  B^)  in  the  disturbing  function  of 
Eq.  (2).  The  unit  vectors  P,  Q,  and  u'  may  be  written  in  terms  of  the 
classical  elements  (Fig.  1)  as 


and 


P 


cos  £2  cos  (j  -  sin  £2  sin  u>  cos  i 
sin  cos  cj  +  cos  i?  sin  cos  i 
sin  i  sin  U 


Q 


-sin  (J  cos  £2  -  cos  <J  sin  £2  cos  i 
-sintJ  sin  £2  +  cos  (J  cos  £2  cos  i 
sin  i  cos  (J 


cos  £2'  cos  u  -  sin,Q'  sin  u  cos  i' 
sin£?'  cos  u  +  cos  £2'  sin  u  cos  i' 
sin  i'  sin  u 


where  u  is  the  argument  of  latitude  of  the  third  body. 


(3) 


(4) 


THIRD  BODY 


V 


Fig,  1  Earth- Centered  Geometry 
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We  ran  express  A  and  B  in  terms  of  the  angles  as 

A  =  P  •  u'  =  cosu(cos  u  cos  AQ  t  cos  i'  sin  u  sin  AQ) 

-  sin  W  cos  i  (cos  u  sin  AQ  -  sin  u  cos  AQ  cos  i ' ) 

+  sin  i  sin  i‘  sin  lj  sin  u 

-sin  <J  (cos  u  cos  AQ  +  sin  u  sin  AQ  cos  i 1 ) 
cos  cj  cos  i  (cos  u  sin  AQ  -  sin  u  cos  AQ  cos  i ' ) 

+  sin  i  sin  i'  cos  u  sin  u  (5) 


B  =  Q  .  u1  = 


whereof?  -Q  ~  Q' . 


In  this  first-order  analysis,  we  assume  that  the  third-body  orbit  is  cir¬ 
cular  because  of  the  small  eccentricities  of  the  orbits  of  the  sun  and 
moon.  The  series  expansions  may  start  from  A  and  B  using  the  com¬ 
puterized  series  expansion  package.  Each  term  of  the  series  may  be 
expressed  in  a  general  form  as 


J  1*  •  J3  *1 


I. 


Term  =  C_  (SI)  (Cl)  (SI3)  (CI3) 

V  • 


sin 

cos 


(JjL  +  J2W  +  J3D) 


(6) 


where 


Cr  T  =  coefficient  of  each  term 
ll*  *  4 

and 


SI  =  sin  i 

L  =  u 

Polynomial 

Cl  =  cos  i 

Angular 

W  =  U) 

Variables 

SI3  =  sin  i' 

Variables 

D  =  AQ = Q 

CI3  =  cos  i 1 

» 

2  2  2  2 

The  two  series  (A  +  B  )  and  (A  -  B  )  in  the  averaged  disturbing  func¬ 
tion  Eq.  (2),  and  other  intermediate  series,  are  expanded  in  literal 
form  and  printed  out  by  a  CDC  7600  computer  at  The  Aerospace  Corpo¬ 
ration  as  shown  in  Appendix  A,  Table  A-l.  The  variations  of  the  classi¬ 
cal  elements  of  the  perturbed  body  (an  artificial  satellite)  may  then  be 
expressed  in  terms  of  the  intermediate  results  shown  in  Appendix  A. 
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da 

dt 


0 


de  15 

jf-  =  -  —  eys  (IMW) 


5|i  =  .  2  -Z£^_L  +  2  e2^  IPO  +  |  e2  (IMO  -  cos  i  IMW) 


d£ 

dt 


3  yr/1  +  32\/cos4 

4  s  l  2  sin  1 

dw  _  3  J 

dt  2  sy[ 


3  ,.2 


5  ,.2 


(A  +  B  )  -  1+2  (A  -  B  )  -  cos  i 


IPS  -  IPC^  +  |  e2  IMS  -  IMcJj 

•]- 


d Q 

dt 


dM 

dt 


=  n  -  2y 


3  2 


!(■  +•■)[} 

'[> 


2  2 
(A^  +  B  ) 


5  , .  2 


-  2  s~y  2  (A  +  B  )  -1+2  (A  -  B  ) 


"H 

] 


5  2  ,  .2  _2. 

4-  e  (A  -  B  ) 


(7) 


where  IMW,  IMO,  etc.,  are  the  intermediate  series  shown  in  Appen¬ 
dix  A,  and 

3 


Rm 


S  =  y/i  .  e2 

n  =  mean  motion  of  the  perturbed  body 
n1  =  mean  motion  of  the  perturbing  body 


Rm  =  mass  ratio,  Rm  =  1  for  solar  perturbation  and 

Rm  =  1/81.3  for  lunar  perturbation 


Those  important  terms  which  give  long-period  and  secular  variations 
of  classical  elements  due  to  third-body  perturbations  may  be  clearly 
seen  from  the  above  equations.  The  averaged  variational  equations  for 
J 2,  J3,  and  can  be  found  in  Ref.  5;  they  are  included  in  Appendix  B 
for  completeness. 

Perturbations  due  to  atmospheric  drag  and  solar  radiation  pressure 
are  not  considered  herein  because  of  their  insignificant  effects  on  very 
long-period  variations  for  higher-altitude  satellites.  The  effect  due  to 
tesseral  harmonics  (J^2  to  ^44’  etc*  )  becomes  important  only  when  the 
mean  motion  of  the  satellite  and  the  earth  rotation  rate  are 
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commensurable.  This  effect  will  be  included  in  Part  II  when  this  first- 
order  theory  is  applied  to  the  prediction  of  the  longitude  of  the  ascend¬ 
ing  node  of  the  nearly  12-hr  orbit  of  GPS  satellites. 

Analytical  Integration 

Those  equations  of  variation  expanded  in  series  shown  in  the  previous 
section  may  be  integrated  analytically  to  yield  a  first-order  solution  if 
the  proper  rates  of  the  angular  variables  are  known.  The  concept  of  an 
intermediate  reference  orbit  or  rotating  ellipse  allows  us  to  have  non¬ 
zero  rates  for  the  Keplerian  angular  elements  Q  and  (J.  The  rates  are 
calculated  from  the  secular  terms  in  the  equations  of  variation.  They 
may  be  summed  up  as 


Q 

=  Q 

+  Q 

+  32  j 

+  £_ 

s 

sun 

moon 

J2 

J 

to 

=  <j 

+  <j 

+  <jT 

+  wT 

s 

sun 

moon 

J 

(8) 


The  secular  terms,  by  definition,  are  those  terms  containing  only  the 
mean  motion,  eccentricity,  and  inclination  which  vary  much  more 
slowly  than  the  node  and  argument  or  perigee.  Those  terms  can  be 
identified  from  the  series,  and  they  are 
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When  the  third  body  is  the  moon,  the  above  equations  should  be  divided 
by  the  earth-moon  mass  ratio  (81.3).  The  above  expressions  agree 
with  those  derived  by  Kozai  .  The  rates  due  to  and  are  the 
secular  part  (terms  without  the  argument  U))  of  Eqs.  (B-l)  and  (B-3). 
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The  mean  motions  n1  of  the  sun  and  moon  are  known,  and  the  integration 
is  performed  on  the  variations  with  secular  terms  removed.  The  inte¬ 
grated  results  for  each  perturbation  (J2,  sun,  etc.)  maybe  expressed 
in  a  series  form  similar  to  Eq.  (6). 


1 


(J1L0  +  J2W0  +  J3Do) 


where 
L 
W 
D 

and 

Lq  -  argument  of  latitude  of  the  perturbing  body  at  t^ 
=  argument  of  perigee  of  the  perturbing  body  at  t^ 
D0  =  the  value  of  ($?-&')  at  t^ 


=  n  (t-t„)  +  Lr 


=  %  <*-V  +  wo 
*  <4  "-‘o'  +D0 


sin 

cos 


(12) 


The  values  of  inclination  and  eccentricity  of  the  perturbed  body  or  the 
satellite  at  the  epoch  time  tQ  are  substituted  into  the  polynomial  vari¬ 
ables  SI,  Cl....  in  Eq.  (12).  The  inclination  of  the  sun,  i  ’ ,  is  nothing 
but  the  inclination  of  the  ecliptic  plane. 

The  inclination  i ' ,  argument  of  latitude  L,  right  ascension  of  the  node 
Q' ,  and  rate  12'  of  the  moon  must  be  computed  in  a  special  way  because 
of  the  earth  equatorial  coordinate  system.  From  Fig.  2  and  the  rela¬ 
tions  of  spherical  trigonometry,  we  have  elements  for  the  moon  re¬ 
ferred  to  the  earth  equatorial  system  as 
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Fig.  3  Variations  of  Inclination  and  Node  of  the  Moon  in 
Earth  Equatorial  Coordinates 


Finally,  the  first-order  solution  is  the  linear  combination  of  the  inte¬ 
grated  variational  equations  due  to  each  perturbation. 
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Comparisons  with  Numerical  Integration 

Orbit  variations  in  terms  of  Keplerian  elements  are  calculated  from  the 
first-order  solution  for  a  span  of  800  days.  The  initial  conditions  at  a 
common  epoch  of  1  July  1985  at  0  hr  are  listed  below. 


0.  01 

5  deg 

a  =  14341.8  nmi,  e  = 

0.  7 

,  i  = 

30  deg 

6  3  deg 

Q  =  40  deg. 


U)  =  60  deg,  M  =  0  deg 
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Results  of  these  six  different  cases  are  compared  with  the  results  com¬ 
puted  by  the  program  ELEMENT,  developed  by  B.  Baxter  of  The 
Aerospace  Corporation^,  which  numerically  integrate  the  singly 
averaged  equations  of  variation  shown  above  in  "Analytical  Integ'  tion." 
The  same  perturbative  model  including  through  J4,  the  sun,  and  the 
moon  is  used  in  the  numerical  integration.  The  comparison  indicates 
that  the  variations  in  eccentricity  and  inclination  predicted  by  the  first- 
order  solution  agree  amazingly  well  with  that  generated  by  ELEMENT 
for  low-eccentricity  orbits  as  shown  by  Figs.  4,  5,  7,  8,  and  9.  The 
discontinuity  in  slope  in  the  curves  of  numerical  integration  is  simply 
caused  by  the  rather  large  step  size  in  plotting.  For  high-eccentricity 
orbits,  the  agreement  is  generally  good  up  to  400  days  as  shown  by 
Figs.  6  through  9.  After  400  days,  the  prediction  by  the  analytical  so¬ 
lution  is  noticeably  different  from  the  true  variation. 


Fig.  4  Eccentricity  Time  History  for  Low  Inclination 
Orbit  (i  =  5  deg) 
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Fig.  9  Inclination  Time  History  for  High  Inclination 
Orbit  (i  =  63  deg) 

This  is  probably  because  of  the  relatively  large  variations  in  the  mean 
orbit  elements  due  to  the  high  eccentricity.  Nevertheless,  the  general 
trend  of  variation  in  the  800 -day  span  has  been  predicted  by  the  first- 
order  solution.  Table  1  shows  the  comparison  for  other  orbit  elements 
Q  and  U  that  are  not  shown  in  Figs.  4  through  9.  Table  2  shows  the 
comparison  for  a  GPS  satellite  with  a  =  14341,  8,  e  =  0.  005,  i  =  45  deg, 
Q  -  265.  4553  deg,  and  (J  =  90  deg. 

The  accuracy  of  the  first-order  solution  may  be  improved  slightly  by 
including  the  eccentricity  of  the  orbits  of  the  sun  and  moon.  The 
series  expansion  for  such  an  inclusion  was  made;  it  is  found  that  even 
to  the  first  order  in  eccentricity  e',  the  total  number  of  terms  in  the 
disturbing  function  as  well  as  in  the  equations  of  variation  has  increased 
nearly  three  times!  The  penalty  in  computing  time  and  complexity 
would  be  too  great  for  a  small  improvement  in  accuracy.  Therefore, 
we  did  not  make  the  comparison  with  numerical  integration. 
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Table  2 


DIFFERENCE  BETWEEN  FIRST-ORDER  THEORY  AND 
NUMERICAL  INTEGRATION 


Days 

e 

x  1  O' 

■  4 

i 

(deg) 

a 

(deg 

) 

U 

(deg) 

f  rom 
Epoch 

N.  I. 

Ae 

N.  I. 

Ai 

N.  I. 

AQ 

N.  I. 

Au 

100 

49.  73 

0. 02 

44.899 

0.  002 

260. *  8 

-0.  01 

91.  95 

-0.  027 

200 

49.  04 

0.  15 

44.827 

0.  002 

255.  55 

-0.  01 

93.  79 

-0. 027 

300 

48.43 

0.  35 

44.  780 

0.  003 

250.  54 

J04 

95.  85 

-0. 022 

400 

47.66 

0.  59 

44.673 

0.  004 

245.63 

-0. 003 

97.  55 

0.  033 

500 

46.  77 

0.  92 

44.661 

0.  004 

240. 65 

0.  03 

99.  56 

0.  109 

600 

45.  90 

1.31 

44.  567 

0.  003 

235. 71 

0.  05 

101. 24 

0.24 

700 

44.  74 

1.  7 

44.  547 

0.004 

230. 78 

00.08 

103.  05 

0.45 

800 

43.  77 

2.  3 

44.480 

0.  001 

225.  81 

0.  10 

104.  71 

0.69 

N.I.  =  numerical  integration  of  the  averaged  equations 
Ae.Ai...  =  (first-order  theory)  -  (numerical  integration) 


PART  II:  APPLICATION 
Background 

The  particular  application,  presented  in  detail  below,  is  for  the  de¬ 
sired  strategy  of  orbit  maintenance.  During  the  mission  design  of  the 
GPS  orbits,  a  family  of  orbits  was  numerically  integrated  with  the 
averaged  equations®  to  study  the  long-term  characteristics.  For  each 
orbit  a  series  of  impulses  or  orbit  corrections  must  be  simulated 
during  the  integration  so  that  the  longitude  of  the  ascending  node  cross¬ 
ing  remains  within  a  specified  tolerance  (±1  deg)  for  repeating  ground 
track  coverage.  The  GPS  satellites  have  a  nearly  12-hr  period,  and 
are  required  to  have  repeating  ground  tracks  every  two  revolutions  at  a 
specified  longitude  of  the  ascending  node  crossing.  The  initial  condi¬ 
tions  are  determined  in  such  a  manner  as  to  synchronize  the  repeating 
ground  track  in  the  presence  of  perturbations.  However,  as  time  goes 
on,  the  perturbations  due  to  the  earth  gravitational  potential  (J2,  ...  J4, 
J 2  2*  •  J4  4).  the  sun,  and  the  moon  drive  the  longitude  of  the  node  away 
from  the  specified  value. 

The  method  of  orbit  maintenance  developed  in  Ref.  12  makes  use  of  the 
scalloping  motion  depicted  in  Fig.  11.  The  perturbation  in  semi¬ 
major  axis  provides  the  restoring  acceleration  for  a  longitude  drift 
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induced  by  a  4v  perturbation.  The  method  functions  open  loop,  how¬ 
ever,  by  utilizing  the  induced  drift  rate  for  the  first  correction  in  all 
subsequent  corrections.  This  leads  to  difficulties  if  the  perturbations 
in  argument  of  perigee  and  right  ascension  of  ascending  node  begin  to 
dominate  the  perturbation  in  semimajor  axis  as  happens  for  certain  of 
the  GPS  orbits.  The  open-loop  technique  also  fails  to  completely  fill 
the  dead  band  even  though  the  perturbed  state,  together  with  the  con¬ 
ditions  of  repeatability,  is  used  to  calculate  each  Av .  This  arises 
because  the  is  not  adjusted  by  observing  the  depth  of  the  scallop. 


Fig.  10  Orbit  Geometry  for  Repeating  Ground  Track 
Coverage  (Q  =  1  orbit) 


An  algorithm  which  applies  the  first-order  theory  for  longitude  predic¬ 
tion,  which  in  turn  provides  the  desired  strategy  for  orbit  maintenance, 
has  been  developed  and  employed  for  simulating  the  orbit  corrections 
during  the  integration  of  the  24  orbits  of  the  GPS  Phase  III  system. 
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A  n  a  I  v  s  i  s 


Let  us  begi  n  wi th  t ho  ea  rlh  equato  r  :a  1  s\sten  .i,  .m-rt.ul  spate.  F,g- 
ure  10  illustrates  the  orbit  geometry.  For  repealing  ground  t  r<o  k 
rovi'rajji's,  we  require  that  the  longitude  lea  rt  )i  - 1  j  \ed  i  of'  the  ast  ending 
node  crossing  of  the  satellite,  AAxX>  repeats  cicry  Q  revolutions,  or 


Qr»,to 


XP>  +  A 


A  XX 


=  Qn, 


XPi 


2X7T 
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I  16  | 


with 

N  =  0,  1.  2 . 

tg  =  epoch  time  when  the  ascending  node  crossing  occurs  at 
longitude  AA  NX 

P  =  nodal  period 

a  =  right  ascension  of  Greenwich 
Q  -  right  ascension  of  the  ascending  node 


Fig.  11  Schematic  Drawing  of  the  Time  History  of 
Orbit  Corrections 
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Assuming  a  constant  earth  rotation  rate  and  a  constant  nodal  regression 
rate  determined  from  the  perturbations  due  to  J2,  etc.,  at  tp,  we  can 
rewrite  the  above  equation  as 


"gV*np“g,xanx  =2V*®'0>NPtf  ,17) 

As  time  goes  on,  the  nodal  period  P  and  the  nodal  regression  rate  Q 
vary  due  to  perturbations;  consequently,  X^NX  sl°wly  drifts  away  from 
the  original  value.  By  varying  Eq.  (17),  we  obtain  the  variation  in 
Aa  NX  as  a  functi°n  °f  the  variations  in  P  and  Q  as 

6XANX  =  ^(V  *  “G]  <5(NP)  +l5^NP  <18) 

Since  6S?  is  not  constant,  the  nodal  variation  should  actually  be 


i0NP 


VNP  .  . 

[£<t)  -fi(t0)]dt 


(19) 


The  term  0(NP)  is  defined  as  the  variation  in  time  measured  in  accumu¬ 
lated  variations  of  nodal  periods  since  the  epoch  tp.  It  may  be  ex¬ 
pressed  as  follows: 


VNP 

<(NP)  =  j  ~  dt 


(20) 


where  6P/P  is  the  variation  in  period  per  one  nodal  period,  which  may 
be  related  to  the  variations  of  the  mean  classical  elements  through  the 
following  derivation. 


Let 


then 


(21) 


6P  6u_ 

P  "  - 

u 


(22) 


where  u  is  the  mean  angular  rate  of  the  argument  of  latitude  in  the  per¬ 
turbed  orbit.  It  has  two  components 


_L  _\ 

u  =  n  +  u 


(23) 
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The  first  is  the  classical  component  which  is  the  mean  motion  of  the 
orbit;  the  second  component  is  the  mean  variation  due  to  perturbations 
(as  denoted  by  \  following  Herrick^).  We  know  from  Ref.  13  (Chapters 
IS  and  17)  that 


_ \  \  • 

u  =  -  Q  cos  i  =  -Q cos  i  (24) 

\  • 

where  Q  or  j?  is  the  mean  (averaged  over  one  period)  nodal  rate  of  the 
perturbed  orbit.  By  varying  Eqs.  (23)  and  (24),  we  obtain  to  the  first 
order  that 


6u  =  in  -  6Q  cos  i 


(25) 


Substituting  Eqs.  (23)  through  (25)  into  Eq.  (22),  we  have 

iP  _  in  -  6  Q  cos  i 
P  n  -  Q  cos  i 


(26) 


•  •  • 

where  6Q  is  simply  Q(t)  -  42(tg),  and  the  variation  in  the  mean  motion  in 
derives  from  the  variation  in  semimajor  axis  ia  as 


in  = 


(27) 


After  combining  Eqs.  (26),  (27),  and  (20),  and  substituting  Eqs.  (19) 
and  (20)  into  Eq.  (18),  we  obtain  the  relation  between  the  longitude 
variation  and  the  variations  in  the  mean  classical  elements. 


iX 


A  NX 


(28) 


where 


QG  ~  ^0* 
n  -  0(tQ)  cos  i 


The  averaged  equations  of  variation  of  S2{t)  are  derived  and  listed  in 
Part  I  of  this  paper.  For  satellites  with  mean  motion  commensurable 
with  the  earth's  rotation  rate,  the  effects  due  to  tesseral  harmonics 
must  be  included.  The  variations  in  semimajor  axis  are  induced  pri¬ 
marily  by  the  tesseral  harmonics.  Therefore,  terms  with  resonant 
angles  must  be  carefully  examined  and  integrated.  Kaula^  provided  a 
set  of  integrated  equations  in  a  general  form  for  perturbations  due  to 
tesseral  harmonics.  For  the  12-hr  (Q  =  2)  orbits,  we  include  J^f  (2,  2,  1, 
1),  (2,2,0,-!)],  J32[(3,2f  1,0),  (3,  2,  2,2),  (3,2,0,  -2)],  J42[  (4.2,  1,  -  1), 
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(4,2,2,  1)],  and  J44[(4,  4,  1,0),  (4,  4,  0, -2)|  as 


¥  -  "£[("/  2(<-2ptql  F^P  c/p,  c/mMJ 

-E  [(?/(%*) 


n 

”T~ 

e  sin 


G/  C  a  (29) 

Zpq  tmpq 


where 


/mpq 


Vc^  +s? 


'/mpq 


Im  4m 


2  2 
C/  +  S; 

Im  X  m 


Aj  =  tan 

/  m 


'&) 


cos  -2p)a»  +  (/-2p+q)M  +  m(i?- a^,)- A^J 
(/- Zp)CJ  +  (/- 2p+q)M  +  m(C-  a^.) 

sin[  (/ -2p)<J  +  (/-2p+q)M  +  m (0-  ] 

(l-Zp)O)  +  (l- 2p+q)M  +  m(fi-aG) 


/-  m  odd 


=  tan 


'(® 


C-m  even 

M  =  mean  anomaly  of  the  satellite 


The  functions  F22j»  ^221’  etc*»  are  t^ie  inclination  and  eccentricity 
functions  tabulated  in  Appendix  C. 

For  orbits  with  periods  of  nearly  12  hr,  the  divisors  in  Eqs.  (29)  and 
(30)  may  become  very  small. 


(l-Zp)O)  +  M  +  Ziff-  <iG)  «  0 


(30) 


Generally  speaking,  the  first-order  solutions  [Eqs.  (29)]  usually  break 
down  in  the  vicinity  of  a  resonance.  Fortunately,  the  orbit  maintenance 
corrections  confine  the  value  of  iCf  to  within  a  small  range  for  repeating 
ground  tracks.  Consequently,  this  stabilizes  the  small  divisors  and 
makes  the  first-order  solution  valid.  In  the  event  the  divisor  becomes 
smaller  than  the  known  angular  rates  M,  <j,  and  G  of  the  orbit 
(»2  X  10*5  rad/day),  the  following  approximation  should  be  adequate 
for  predicting  the  variations  from  about  200  to  400  days  for  orbit  main¬ 
tenance  purposes. 


constant 


constant 


(31) 
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The  right-hand  side  of  the  above  equations  can  be  obtained  by  differen¬ 
tiating  Eqs,  (29).  The  additional  nodal  rates  from  Eq.  (31)  should  be 
added  to  the  total  nodal  rate  &s  as  shown  in  Eq.  (8). 

An  Algorithm  for  Orbit  Maintenance 

The  desired  strategy  of  orbit  maintenance  is  to  keep  to  a  minimum  the 
number  of  corrections,  i.e.,  impulses,  within  a  given  time  span.  In 
other  words,  the  4V  should  be  correctly  determined  so  that  the  time 
interval  between  two  adjacent  corrections  is  a  maximum  without  getting 
outside  the  specified  range  as  shown  by  Fig.  10.  For  GPS  satellites, 
the  time  interval  between  corrections  varies  from  200  to  400  days;  it  is 
too  costly  to  numerically  integrate  the  orbit  from  tj  to  t2  (Fig.  11)  in  an 
iterative  manner  for  the  determination  of  4V  at  tj.  Thus,  the  first- 
order  solution  which  predicts  the  longitude  variation  due  to  the  effects 
of  the  sun,  moon,  and  earth  gravitational  potential  as  shown  by  Eq.  (28) 
is  employed  to  determine  the  Av  at  each  orbit  correction. 

After  the  integrations  in  Eq.  (28)  are  performed,  the  longitude  varia¬ 
tion  may  be  expressed  in  series  form 


N 

- 

sin 

sin 

- 

dAAIVTV  =7  b. 

A  NX  1 

i=l 

cos 

[yj(t-t.)  +g£]  - 

cos 

(gi> 

where  ^ANX  *s  the  variation  of  longitude  of  ascending  node  from  the 
epoch  value  at  time  tj  when  an  orbit  correction  is  made,  and  B^,  >^,  and 
are  constants  computed  from  mean  orbit  elements  i,  e,  n,  and 
initial  conditions. 

As  shown  by  Fig.  10,  the  first  orbit  correction  should  be  made  at  tj 
when  the  value  of  ^ANX  reaches  the  lower  bound.  The  desired  strategy 
requires  that  the  maximum  variation  of  ^anX  occurs  at  tm  (the  re¬ 
storing  force  is  the  tesseral  harmonics)  with  its  magnitude  nearly  equal 
to  the  specified  tolerance,  say  2  deg.  Let  tj  be  the  time  of  an  orbit 
correction  and  tm  at  some  time  between  tj  and  tj+|,  which  is  the 
time  of  the  next  correction.  We  have 


N 


iXANX(tm! 


=  2  deg 


-E 

i=l 


B. 


sin 


*il  ' 


(gt) 


cos, 


cos 


With  the  slope  equal  to  zero  at  t  ,  we  require 


A  NX 


N  / cos  \ 

(t  )  =  B.y.  (  |[y. (t  -t.)  +  g. 

m'  /  -■  11  l  /  1  m  y  1 

i  =  l  IVsin/ 


]  +  4A.  =  0  (34) 


The  second  term  on  the  right-hand  side  of  the  above  equations  is  intro¬ 
duced  when  an  orbit  correction  is  made  at  tj.  The  introduced  constant 
rate  AX  •  at  t-  is  an  unknown  parameter  to  be  determined.  With  the 
above  .equations  we  should  be  able  to  determine  the  two  unknowns  t 
and  4A-.  An  interative  formula  derived  from  Newton's  iteration 
methoa  is  constructed  to  compute  t  . 


(tmo) 


~  tmA  ~  F'7t 
0  /  m 


where 


N  f  f  sin 


F(tm>  =EBi 


i=l  LIcos 


g.  t  y.  (t  -t.) 

Bi  'i'  m  y 


di  "X<,Xbound 


=Evf(tm-V 


6.  =  y.  (t  -t.)  +  g. 
i  i  m  j  ei 

^Abound  =  sP€Ci^e(^  bound,  2  deg 

The  parameter  X  is  a  factor  which  is  slightly  smaller  tjian  1,  so  that 
the  true  variation  will  not  exceed,  the  specified  bound  $At,oun(j.  Once 
t  is  obtained  after  iteration,  A\:  is  determined  from  Eq.  (34).  By 
differentiating  Eq.  (28),  we  can  determine  the  correction  for  semi¬ 
major  axis  from  the  value  of  ^Aj  as 


(¥),  - ! 
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or  a  corresponding  correction  to  the  period  of 


(38) 


The  initial  guess  of  t,^  may  be  determined  from  an  approximation 
the  small  divisors  in  ®Eqs.  (29)  or  (30)  vanish  and  6Q  is  zero.  The 


Usually,  the  longitude  of  the  ascending  node  crossing  is  sustained  to  the 
initial  value  within  a  dead  band  of  ±  1  deg.  The  sustenance  strategy 
functions  by  applying  impulses  at  perigee  when  the  longitude 
violates  a  boundary  of  the  dead  band.  The  impulse  is  performed  to 
adjust  the  period  in  order  to  give  a  repeating  ground  track  at  the  time  of 
the  correction;  that  is,  the  synchronous  period  is  redetermined  to 
include  the  effect  of  perturbations  in  the  elements  and  the  predicted 
variation,  using  the  new  algorithm  described  above.  The  new  semi¬ 
major  axis  is  determined  by  the  following  relation: 


*  4p,  ♦  Ap2  +  4p3  =  Prequired  (■ 

where 

a  =  new  semimajor  axis 

=  period  variation  predicted  for  optimum  orbit  sustenance 
(the  new  algorithm) 

^required  “  27y4*G  " 
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and 


Of,  =  rate  of  the  right  ascension  of  Greenwich, 

°  7. 292 11 543 IX 10 "5  rad/sec 

• 

Q  -  mean  rate  of  the  right  ascension  of  the  ascending  node 

The  period  variation  due  to  J?  effect,  AP?,  is  obtained  from  Claus  and 
Lubowe15. 


C  •  2  . 

-5  sin  1 


(1+e  co  scj ) 


(42) 


with  a»  =  earth  equatorial  radius,  r  =  a(l-eZ)/(l+e  cosu) ,  and  P  = 
Zn/yJ/HJaX .  The  period  variation  due  to  the  sun-moon  perturbation  is 
similarly  derived  as 


AP-, 


9/2 

e“)  Rm  cot  i  (xy1  U1  +  yy'  U2) 


(43) 


where 


Rm  =  mass  ratio  = 


l:solar  perturbation 
1/81.  3;  lunar  perturbation 
.  2  i' 


=  cos  (L.-AS?)  -  2  sin^  ^  sin  L  sin  AS? 


=  (i-2  sin2  jj 


l[sin(L-4X?) -2  sin2  ^  sin  L  cos  AS? ] 

+  sin  i  sin  i*  sin  L 

y1  =  dy/di  =  cos  i  sin  i'  sin  L  -  2  sin  i  sin  (L -ASh 
+  4  sin  i  sin2  ^  8*n  L1  COB  AS? 

ZTt 

=  J"  sin  u  cos  u  du/[l+e  cos(u-W)]^ 

0 

■  2 lTTe2  sin2w(l  +  e2/2)/[4(l  -  e2)11/2l 
2  7T 

/2  6 

sin  u  du/[l+e  cos(u-u)] 

0 

=  [>r/(l-e2)^^2][l  +  (5-^pcos2w)e2  + -^  (5-7 cos2o»)e^] 


U1 


U2 
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For  GPS  satellites,  the  orbits  are  nearly  circular  (e  ss  0),  and  Eq.  (43) 
can  be  simplified  to  the  following  form: 


“p1  =0)(^~)  (^)  Rm  COt  1  (Yy,)  (44> 

This  algorithm  has  been  incorporated  into  the  computer  program 
ELEMENT  for  an  analysis  of  the  orbit  perturbations  for  the  Phase  III 
GPS  satellites. 

Figure  12  gives  an  example  of  orbit  maintenance  with  this  algorithm. 
This  example  is  one  of  the  12  cases  at  55-deg  inclination  of  the  Phase  III 
system^.  For  this  particular  satellite,  each  correction  requires 
0.6  ft/sec  or  less  in  4v,  and  the  total  Av  used  in  the  10-yr  period  is 
about  7.  5  ft/sec.  The  first  and  third  corrections  are  caused  by  a  slight 
overshoot  at  the  upper  bound  (271  deg).  This  occasional  overprediction 
seems  to  reveal  a  limitation  of  the  first-order  theory  assumed  in  the 
new  algorithm.  To  minimize  the  number  of  over-  and  under-predic¬ 
tions,  a  self-adjusting  mechanism  is  incorporated.  After  an  overshoot 
correction,  as  shown  by  1 3  in  Fig.  12,  the  predicted  maximum  value  of 
the  variation  in  in  the  next  interval  is  scaled  down  by  10%.  If  the 

previous  maximum  oX^NX  *s  rnore  than  10%  undershoot,  the  next  maxi¬ 
mum  MaNX  *s  sca^ed  up  proportionally.  This  self-adjusting  scheme  is 
believed  to  be  more  efficient  in  computation  than  actual  inclusion  of 
those  second-order  terms  of  the  perturbations  which  would  significantly 
slow  down  the  computation. 

The  long-term  (10-yr)  variations  of  the  classical  elements  of  the  12 
GPS  satellites  and  the  orbit  maintenance  history  may  be  found  in  Ref.  16. 
Results  of  those  variations  revealed  the  interesting  fact  that  the  value 
of  inclination  increases  when  the  node  Q  is  between  0  and  180  deg,  and 
the  inclination  decreases  when  the  node  is  between  180  and  360  deg. 

The  drift  rate  of  inclination  vanishes  when  the  node  is  around  0  and 
180  deg.  With  the  analytical  expansions  derived  in  Part  I,  we  are  able 
to  see  that  for  this  type  of  orbit  (e  as  0)  the  perturbations  on  inclination 
are  mainly  due  to  the  sun  and  the  moon.  Neglecting  terms  of  e^  or 
higher,  we  obtain  the  averaged  variational  equation  in  inclination  due  to 
third-body  perturbations  from  Eq.  (7)  as 


in(i?  -  £3) 

(45) 

Those  parameters  with  a  subscript  3  are  of  the  third  body.  The  node 
of  the  sun  is  zero  and  the  node  of  the  moon  is  confined  to  within  ±15  deg 
on  the  earth  equatorial  plane  where  Eq.  (45)  is  derived.  With  this 
equation,  one  can  explain  why  the  variation  of  inclination  is  a  function  of 
the  right  ascension  of  the  ascending  node  Q  . 


2  3 

di_3^/^3\n 
dt  4  n  \ r3  / 


cos  i  sin  2  i^  s 


+  2  sin  i  sin^  i^  sin  2 (Q - 
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Fig.  12  Variations  of 


under  Orbit  Sustenance 


The  above  example  clearly  indicates  that  the  analytical  solution  has 
provided  insight  into  the  long-term  variations  of  the  orbit  elements  of 
the  GPS  satellites.  This  is  particularly  important  for  mission  design¬ 
ers  to  understand  the  long-period  characteristics  of  the  orbit.  Figure  13 
shows  how  closely  the  first-order  solution  predicts  the  inclination 
variations  of  three  GPS  satellites  in  three  different  orbit  planes.  Even 
after  10  years,  the  error  is  only  slightly  over  10%.  The  computer  time 
required  to  generate  the  30  points  of  prediction  is  about  15  sec  (CDC 
7600). 
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p 


Fig.  13  Inclination  Variation  of  the  Phase  HI  GPS  Satellites 


CONCLUSIONS 

The  equations  of  variation  perturbed  by  the  sun -moon  attractions  are 
expanded  by  computer  in  terms  of  the  familiar  classical  elements.  The 
motions  of  the  sun  and  the  moon  are  represented  by  circular  orbits  with 
the  nodal  rate  of  the  moon  approximated  by  the  projected  mean  value  in 
the  earth  equatorial  system.  Results  indicate  that  the  first-order  solu¬ 
tion,  analytically  integrated  with  the  secular  part  included  in  the  refer¬ 
ence  orbit,  can  successfully  handle  the  coupling  between  the  sun-moon 
attraction  and  the  oblateness  effect.  Good  accuracy  in  predicting  the 
variations  of  orbit  elements  up  to  800  days  has  been  demonstrated  with 
examples  of  the  GPS  satellites.  However,  such  accuracy  of  orbit  pre¬ 
diction  is  limited  to  100  days  for  high  eccentricity  orbits  (e  =  0.7). 

This  is  primarily  because  of  the  relatively  large  perturbations  due  to 
J2  at  high  eccentricity  which  make  the  second-order  terms  important. 
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An  algorithm  applying  the  first-order  solution  has  been  developed  to 
achieve  the  desired  strategy  of  orbit  maintenance  for  the  GPS  Phase  III 
system.  The  result  of  a  10-yr  simulation  shows  that  the  total  Av  re¬ 
quired  to  maintain  the  longitude  within  *  1  deg  is  around  7.  5  ft/ sec  with 
each  correction  of  about  0.  5  ft/sec.  This  algorithm  will  be  useful  in 
determining  the  time  and  magnitude  of  the  orbit  correction  maneuvers 
for  the  GPS  satellites.  The  analytical  solution  has  helped  to  clarify  the 
long-term  variations  of  the  GPS  satellite  orbits. 
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APPENDIX  A 


The  package  of  series  expansions  used  in  this  analysis  is  the  latest 
version  of  the  computerized  series  expansion  package  developed  by 
R.  Broucke  of  the  University  of  Texas.  After  the  terms  are  properly 
combined  by  hand,  the  two  series  may  be  rewritten  as  shown  in  Table 
A-2.  This  compact  series  may  be  useful  for  other  types  of  analysis 
such  as  deriving  the  equations  of  variation  from  the  averaged  disturbing 
function  [Eq.  (2)]  in  the  Delaunay  variables.  However,  the  results  of 
Table  A-l  will  be  used  for  continued  series  manipulation  by  computer 
to  avoid  any  possible  human  error.  We  can  then  derive  those  partial 
derivatives  of  (A^  +  B^)  and  (A^  -  B?-)  with  respect  to  the  classical  ele¬ 
ments.  Again,  the  operations  were  performed  by  the  computer;  results 
are  shown  in  Tables  A-3  through  A-5.  The  nonzero  series  are  denoted 
by  the  following  symbols: 


IMW 

IPS 

IMS 


_  d(A2  -  B2) 

=  a(AZ  4  Bl) 
a  sin  i 

=  a(A2 .  b2) 

d  sin  i 


IPO 

I  PC 

IMC 


-  a(A2  +  B2) 
dQ 


=  a(A2  4  B2) 
acos  i 


a(A2  -  b2) 

9cos  i 


IMO 


-  a(A2  -  B2) 
dQ 


(A-l) 


Aerospace  internal  correspondence.  Not  available  for  external 
distribution. 
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Table  A-l 


COMPUTER-GENERATED  SERIES  FOR  (A2  +  B2)  AND  (A2  -  B2) 


♦(♦i-l/2*S13«*i-l/2*SI**2*3/**Sl**2*S13**2  3 
♦<*1*SX*C1*SX3*C133*C0S<  0  3 

♦  <♦ 1/**S I**2*S I  3**2 3 *CUi (2*03 

♦  (♦iM*SI**2U/**SI**2*.; I3-i/d*SI**2*SI3**23*COS  <2*1-2*03 
♦<-l/2*SI*Cl*Sl3-i/2*Sl*Ci*S13*C133*CQS C2*L-0( 

♦(♦l/2*SI3**2-3/**S I**i*  S13**23*CQS(2*L 3 

♦  <♦*/ 2*SI*CX*S 1 3-*/ 2* il* C 1*51 3*C 1 3  3*CQS (2*L*0I 

♦  (♦i/4*Sl*«2-i/ **51**2*,: i3-l/«*Sl**2*SI3**23*CQS<2*l*2*03 

*2-S2> 

♦(♦l/**S13*«2-l/**Cl*SlJ**2-l/a*Sl*«2*S13**23*CQS(2*M-2*D3 

♦  <♦ 1/ 2*5 1*5  1 3*C 13-1/2*51 *tI*S13*C13)*C0S(2*M-03 
♦<*l/2*5l**2-3/**Sl**2*i 13**2 3*CCS  <2*U X 
*<“1/2*51*5 I3*CI 3-l/2*Sl *C  1*S  1 3*C  I  33  *COS  I  2*m*03 

♦  <*1/M5I3**2*1/**C1*SIJ*«2-1/8*SI**2*S13**2  3*C0$<2*U*2*0) 

♦(♦l/**l/**C13-l/8*Sl3«*2*X/**Cl*l/**CI*CI3-l/a*Cl*S13**2-l/8*Sl**2-X/a*Sl**2*CI 

3*l/16*SI**2*S13«*2l*C05(2*L-2*h-2*OS 

♦(♦1/**SI*5 I3*1/**S1*SI3*C13*1/**S l*C 1*S13* 1/4*5 I*C 1*5 I3*C 13 >*C0S< ?»l -2*«-D3 
♦(♦3/8*Sl**2*S13**2)*CUS  <2*1-2** 3 

♦  (♦l/4*$I*$l3-l/**SI*5i3*tl3-lM*5X*Cl*513*l/**SI*Cl*S  I3*tI3)*CO$(2*L-2*W*D) 

♦<♦1/4-1 /4*CI 3-1 /«*31J**2-l/4*Cl*l/4*Cl*Ci3*l/8«CI*S13**2-l/8*$I**2*l/e*SI**2*Cl 
3*1/ 16*S l**2*S 13** 23* CCS  <2* L-2* **2*0 1 

♦(♦l/4^1/4*C13-l/t»*S 13**2-l/4*Cl-l/4*Ci*Ci3^1/8*CI*SI3**2-l/8*SI**2-l/2*SI**2*CI 
3*l/16*SI**2*SI3**23*C0$<2*L*2*W-2*03 

♦<-l/4*SI*SI3-l/4*SI*513*Cl3*i/4*SI*Cl*Sl3^1/4*Sl*Cl*5l3*Cl3)*C05(2*l^2*tt-0> 
♦(♦3/8*SI**2*S13**23*Cu$ <2*l*2*H3 

♦<-i/4«$I*S I3*1/4*SI*513*CI3-1 /**SX*C1*SI1»X/**SI*CI*SI3*CIS3*C0S<2*1*2***03 

♦<n/*-x/4*ci3-i/**>i3*.2*i/**ci-i/a*ci*cx3-i/a*ci*si3**2-i/a*si**2*x/a*ji**2*ci 

3*X/16*SI**2*SX3**23*C0S<2*L*2*«*2*03 

A2  +  B2  =  A2  +  B2 
A2  -  B2  =  A2  -  B2 
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Table  A -2 


2  +  B2 


SERIES  OF  (A2  +  B2)  AND  (A2  -  B2) 


1  2  2  12  2 
^  (1  +  cos  i)  (1  +  cos  i^)  +  2  s*n  *  s*n  *3 

1  2  2 

+  Tr  [sin2  i  sin2  cos45  +  sin  i  sin  cos2  45] 

L  1  /.  3  .  2  A  .  2  .  , 

+  •^■11-^  sin  if  sin  i^  cos2  u 

+  ^  sin2  i  [(1  +  cos  i)2  cos  2(u  -  AQ) 


2 

+  (1  -  cos  i)  cos  2  (u  +  45)] 

+  ^  sin  2  i  sin  i^  [(1  -  cos  i^)  cos  (2u  +45) 

-  (1  +  cos  i,)  cos  (2u  -45)] 

1  2/3  2  \ 

2  sin  i  II  -  2  s^n  *3  J  cos  ^  (*> 

+  ^  sin2  ij  sin  i  [(1  -  cos  i)  cos(2c»;  -45) 

-  (1  +  cos  i)  cos  (2W+45)] 

+  ^  sin2  i^  [(1  -  cos  i)2  cos  2(u> -AQ) 

+  (1  +  cos  i)2  cos  2  (<J  +  45)] 

3 
8 


2-B2 


3  2  2 

+  a  sin  i  sin  i,  [cos  2(u-a>)  +  cos  2  (u+oi)] 


+  ^  sin^  [sin^  j  cos  2  (u-W+45)  +  cos*^-  cos  2(u+w+45)] 


4  i 


1  4  A  * 

+  ^  cos  -^-[sin  ^-cos  2  (u+w- 45)  +  cos**i  cos2(u- a» -45)] 
+  sin  i  sin  i^  sin2  [sin2y  cos  (2u  -  2W+45) 

-  cos^|  cos  (2u+2w+45)] 

+  sin  i  sin  i,  cos2  -y-  [cos2  y  cos  (2u  -  Zu-AQ) 


4  i 


-  sin2 y  cos  (2u  +  2i»>-45)] 
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Table  A-3 


SERIES  OF  IPO  =  d(AZ  +  B Z)/dQ,  IMO  =  d(A2  -  BZ)/dQ 
AND  IMS  =  d(A2  -  BZ)/d  sin  i 


IfD* 

♦<-X*Sl*CI*SA3*C13l*SlN(  U» 

*<-l/2*SI**2*SX3**2X*SIi«  <2*01 

♦^♦l/2*SI♦•^♦l/2♦Sl**2•: I3-1/**SI**2*SI3**2»*SINU*L-2*DI 

♦  l-l/2*SX*CI*SX3-X/2*Sl*Cl»S13*CI3l*SIM2*L-i>> 

♦<-i/2*$I*C I*S13*1/2*SI*C I*SI 3*Cl3 >*SIM  2*L*0 1 
♦(-l/2*SI**2*i/2*$l**2*C X3*iM*Sl**A*SX3**2l*SlNX2*t.»2*D) 

mo- 

♦^♦l/2♦SI3*•2-l/2♦CI•Sli**2-W*•SI••2*SI3**2 l*SINC2**-2*D) 

♦  (  *1/ 2*  Sl*S 1 3  +  C 1 3-1/ 2 *S 1 *b i*S 1 i*C l 3l*i iN (2*K-D) 
♦(♦l/2*SI*SI3*CX3*l/2*iI*CI*SI3*CI3l*SIN(2*K*0> 

♦l-i/2*SI3**2-l/2*CI»SI3  ♦*2'*1/4*SI**2*SI3**2»*SINX2*k*2»I)I 

♦<♦l/2♦i/2•Ci3-.^/**$l^*•2♦i/2♦Ci♦l/^*CI*CIi-^/*•CI•SIi••2-l/*•SI*♦2-l/*•SI••2•Cl 
3».i/e*SI**2*SI3*»2J*i IM2*L-2*W-2*0> 

»««*/^*SI*SI3*A/^*Si*ii3*ti3*i./4*Sl*Ct*SI3*l/MSI*CI*$I3*CI3»*SINX2*L-2*W-D» 

♦  (-AM*SI*SI3»l/**SX*SiJ*CX3*i/<>«I*Ci*Si3-l/**SI*Cl*Sl3*Cl3»*SlNI>*».-2*W*D> 

♦  <-i/2»X/2*Ci3*i/MSIi**i*i/2*Cl-i/i*Ci*CI3-i/MCI*SX3**2*i/4*SI**2-l/«*SlM2*CI 

3-X/e*$I*«2*513*»2l*S  XM2*L-2«*«2*i>> 

♦(♦l/2^1/2*Ci3~X/**SI3**2”l/2*CI”i/2*CI*Ct3*i/%*CI*SI3**2-X/^»SI**2-l/**SI**2*CI 
3»i/8*SI**2*SI3**2)»i lM2»l*2*M-2«n) 

♦  «-l/**Si*SI3-XM*SI*SIl*U3*i/**SA*Cl*SII*i#**SI*Ci*S  X3*CI3I*SIN«  2U«2*W-DI 
H«l/4*SI«SI3-l/4*$l*SiJ*CX3*Xf4*SX*CI*SI3-X/MSI*CI*SI3»Cl3>»StN<2*l»2*W»D> 
♦<-i/2*i/2*C13*X/**SI3**2-l/2*CI*l/2*Cl*CI3»l#**CI*SI3*»2*l/**SI**2-lM*SI**2*CI 

3-l/«*SI**2*Sl3**2l*SIM2*l*2*<M2«-0» 


IHS- 

♦  (♦X/2*SI3*CX3-i/2*Ci*Sl 3*C 13 » *C0S  I2*«-D> 

♦  l-i*Sl*SI3**2)*CM  12**1 

♦|-1/2»SI3*CI3-1/2*CX*SI1*CX3I*CQS12***D> 

♦1*X/** JI 3*l/**iX3*C  X3*X /**Cl*S 1 3*  1/**CI*S 13*C X3X*CQS ( 2*t-2*M-0l 

♦  l»l/2*Sl*Sl3**2l*C0SU*L-2**> 

♦  <*X/<.*SIi-l/**S13*CI3-i/<.*U*SI3*l/**CI*SI3*Ci3l*CUSC2*l-2*t»*0> 
♦<-l/^$U-l/»*X13*CII*l/^CX*SJ3*lM*Cl*SU»Cl3»«COJU*l*X*»*OI 
♦(♦i/2*3I*Sl3**23*CUil2*l*2*KI 

♦  4-1 /**  S  13*1  M*S  1 3*C  1 3-i /**CI*SI3*X/**CI*SX3*CI3I*C0S( 2*C*2*W*  0) 
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Table  A-4 


SERIES  OF  IPS  -  0<A2  +  B2)/d  sin  i,  IPC  -  d(A2  +  B 2)/d  cos  i 
AND  IMW  -  d(A2  -  B 2)/dui 


IM- 

♦(♦X*SI*S13**2! 

♦(♦x*ci*si3*ci3i*cOi(oi 

♦  4-i/2*Cl*SI3-X/2*CI*3U*CI31*CQS  1 2*1-01 

♦  4-l*SX*Sl3**2l*COSU*U 

♦(♦I/2*CI*SI«-l/2*Cl*>lJ«C13l*C03(2*t*0» 

l*C- 

♦(♦i*Cl-l/2*CI*S13**2» 

«(*X*3I*S(3*C(3>*C0$(DI 
♦l-l/2*CI«S I3**2l*COS42*01 

♦  l-i/2*Cl-l/2*CI*Cli*lM*Ci*SI3**21*C0S  12*1-2*0) 

♦  4-l/2«Sl*Si3-l/2*SX*SU*C13)*C0S<  2*1-0) 

♦<*l/2«CI*SI3**2l*C0S (2*11 

♦ C*i/2*SI*S 13-1/ 2*i 1*511 *C1 31* COS (2*1*0! 

♦l-X/2«CI*X/2*CI*C13*i/**Cl*$13**2  )«C0S 4 2*L*2*0 1 

IHk* 

•(-i/2*S13**2*X/2*Cl*3lj**2«l/4*Sl«*2*S13**2>*SXN(2*M-2*i>) 

♦  <-l*Jl*SI3*Ci3*l*SI*CI*SI3*CI3)*$(N(2*t*-0) 

♦("i*Sl**2*322*Sl**2*Sl3**2)*SIN(2*a) 

*(*X*3I*$I3*CI3«X*SI*U*$13*CI3)*S1N<2*«*0) 

*<-l/2«SI3**2-l/2*Cl*S13**2*X/4*Sl**2*S13**2t*SlN(2»U*2*D) 

♦  <*l/2*l/2*Cl3-X/**SI3**2*l/2*Ci*l/2*Cl*Cli-lM*CI*Sl3**2-XM*SI**2-lM*$I**2*Cl 

3*X/S*SI**2*Sl3**2)*SlN(2*L-2*8-Z*0) 

♦  ( *l/2*Sl*S 14*1 /2*S X*SIi *013*1 /2*SI*C1*S13*X/2*SI*CI*SI3*C13)*SINI 2*1-2* h-01 

«( *3/4*SI**2*SI3**2 )*SX« ( 2*1-2** I 

♦(♦i/2*SI*SI3-l/2*Sl*S|j*C13-l/2*Sl*CI*S13*X/2*SI*CI*Sl3*CI3)*SIH(2*l-2***D) 

♦(♦l/2-l/2*Cl3-l/**SI3*»2-l/2*CI*l/2*C  l*C  1 1*1M*CI*$  13**2-X/4*SI**2*X/4*SI**2*C  I 
3*ISa*SI«*2*Si3**2)*SXN(2*l-2*W*2*0) 

♦(-X/2-l/2*CI3*i/4*SIi**2*i/2*CX*X/2*CI*CI3-l/«*Ci*SI3**2*l/4*SI**2*l/**$X**2*Cl 

3-X/8*SI**2*S  X3**2 )*S 1N(2*L«2*W-2*D 1 

♦  (♦X/2*SI*SI3*l/2*$I*S13*C13-i/2*Si*CI*SI3-X/2*$l*CI*$I3*CX3»*SIN(2*l*2*(|-|i) 

*(>3/4*SI«*2*SI3**2)*S1’M2U*2*m) 

*(*X/2*SI*Sl4-X/2*SI*SiJ*Cl3*X/2*Sl*CI*Sl 3-X/2**l*CX*S I3*CI3)*SI«M 2*L*2*M*0) 
*<-l/2«X/2*CI3*X/4*S13**2-X/2*CI*X/2*Cl*Ct 3*X/«*CI*X13**2*X/4*SI**2-X/4*SI**2*CI 

3-l/8.*Sl*«2*S13**2  l*S  1M2*L*2*U«2*DI 
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Table  A-5 


SERIES  OF  IMC  =  d(A2  -  B 2)/d  cos  i 


ISC* 

*(-l/**SI3**2*l/<t*Cl*»iJ**2  >*C0S(2*¥-2*D» 

*(-X/2*SI*SI3*CI3>*C0S(2*M-(» 

*(-l*CX*X/2*Cl*il3**2 

*(-l/2*SX*$l3*C131*C0i(2*ii*0) 

»(*X/**SX3**2*l/4*Cl*SI3**2l*C0S(2*¥*2*0> 

♦  t  *XM*i  M*C  U-* /#*S  I3**2*l  /**Cl*iM*Cl*CI  3-1  /tf*Cl*S 13**2 >*COS  (2*1-2  *¥-2*01 
»(*l/**SI*Sl3*i/4*Si*S13*C13>*C0S(2*L-2*M-0) 

♦<-i/4*C l*$ 1 3**  2 1  *C  OS  (2*1-2  •¥  I 
*(-XM*SI*S l3*4/**Si*Sl3*CI3l*C0S(2*L-2*¥*0l 

M-X/**l/**CI3*i/e*S  13*«2*l/**Cl-X /**CI*C13-1/9*CI*S1»**2)*C0$  (2*1-2* W*2*0> 
*(-X/*-X/**Ci3*X/0*S13**2*XM*Cl*X/4*Ci*C13-X/O*Cl*Si3**2l*C0S (2*1*2 *¥-2*0 1 
♦(♦X/4*SI*SI3*l/**Si*SiJ*CI3)*CQS( 2*t*2*¥-0> 

*(-l/**CI*SXi**2»*C0S (2*  l*2*¥ I 
*(-lM*Sl*SI3*i/**Sl*SU*C13)*CU$(2*l*2*¥*Dl 

*(*1/4-X/**CX3-a/#*S 13**  2*1 /4*Ci-* /**C1*CI 3-l/¥*CX*S I3**2 )*C0S (2*L *2*¥*2*0 > 


Note  that  the  designed  package  cannot  take  partial  differentiation  direct¬ 
ly  with  respect  to  the  inclination  i,  which  is  treated  as  a  polynomial 
variable.  Thus,  the  intermediate  parameters  sin  i  and  cos  i  are  used 
with  the  following  relations  to  complete  the  differentiation. 

dF  dF  dF  .  .  , . 

rrr—  =  -35 — : r-  COS  1  -  -m - r  Sin  1  (A-2) 

Cl  csin  1  CCOS  l  '  ' 

Note  also  that  the  printed  series  are  in  FORTRAN  language;  they  can  be 
punched  out  and  directly  inserted  into  FORTRAN  programs. 

APPENDIX  B 

J2  Perturbation 


de  45  2  2.  .  2  .  .  _ 

jj-  -  -  jj  n€  2  e(l-e  )  sin  1  sin  2 u) 


di  _  45  2  2  .  ,  . 
dt  "  R£2  n  *  8ln  2  1 


‘( 


14  .  2  . 

TJ  -  sin  1 


da  _  3  9  ,2 

jj-  =  -  2"  n  «2  cos  1  ~  4  n  *  2  008  1 


(||  -  sin2i) 

j  sin  2u) 

[3  5  .  2  .  / 1  5  .  2  .\ 

[2  -5-sln  l+(5-  +  24  81n  f 


( 7  5  .  2  A  2  ‘ 

Iy2  "  4  81n  1je  cos  2U> 


(B-l) 
(cont. ) 


THIS  PAGE  IS  BEST  QUALITY  PEACH CAftLI 
t&M  COPY  FIMUSHJSD  10  DDC 


3?  =  f  n€2  (2  -  I  Sin2  •  In€2  I  2  •  T-sinZ  1  +ll8in4  1 

+  V  (7  ‘  I  si"2  1  +  T"  si°4  +  COS4 1<f>’  [(7  '  'T8in2  si 

+  e2  ^7  +  5  sin2  i  -  ^sin4  i^J  | 


.  2  . 
sin  l 


where 


*2  =  J2  (Re/P2) 

p  =  a(l-e2) 


=  equatorial  radius  of  earth 
J3  Perturbation 

1  2  2 

—  -  n€0  (1-e  )  cos  CJ  sin  i  (1-5  cos  i) 
at  3 

,  *  2 

=  n63  e  coswcos  i  (5  cos  i  -  1) 

jq  2 

=  n«3  e  sinucot  i  (15  cos  i  -  11) 

=  n«3  ^  +e4e  —  sin  (a)  sin  i  (5  cos2  i  -  1)  -  ^cos  i 
where 

<3  =  3/8  (Re/P)3  J3 
J4  Perturbation 

^  =  -n«4  e(l-e2)  sin  2wsin2  i(7  cos2  i  -  1) 

=  i  n€A  e2  sin  2(Jsin  2  i  (7  cos2  i  -  1) 
at  c 

^  =  n«4  cos  i  |2  (7  cos2  i  -  3)  +  e2  [7  cos2  i  -  1 
+  4  sin2W(7  cos2  i  -  4)]  | 
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L 


(B-l) 


(B-2) 


(B-3 ) 
(cont.  ) 


dCJ  ,  „ 
dt  4 


O  A  “7  ?  o 

8-28  sin'1  i  +  21  sin  i  -  sin  w  sin  i  (7  cos^  i  -  1) 


2  I  ,  .  2  .  ,63  .  4  .  ,  .6 

+  e  ^  6  -  14  sin  i  +  -g-  sin  i  +  sin 


U)(b  -  35  sin2  i  +  ^sin*  i)j| 


(B-3) 


where 


t,  =  15/32  (R  /Pr  J, 


APPENDIX  C 


Table  C-l 


221 


321 


322 


311 


330 


421 


440 


211 


THE  F 


I  AINU  lii 
imp  Xpq 

Isin2!,  F220  =  |  (1+cos  i)2 


AND  Gi  _  FUNCTIONS 


15  2 

-g-sin  i  (1-2  cos  i  -  3  cos  i) 


15  2 

-  -g-sin  i  (1+2  cos  i  -  3  cos  i) 


sin2  i  (1+3  cos  i)  -  ~  (1+cos  i) 

15  /,x  m3 
-g-  (1+cos  l) 

105  .  2  .  .  ...  ..  15  ...  .  .2 

— g-  sin  l  cos  i  (1+cos  i)  — g-  (1+cos  i) 

105  ..  .  .,4 

-Jg-  (1+cos  l) 

1. 5e  +  1. 6875e3  +  2. 0390625e5  +  2. 3289388e7  +  2. 587323e9 
+  2.822341eU  +  3. 03933758e13  +  3.2418959e15  +  3. 43255448e17 
+  3.6131875e19  . 


G20-  1  =  ■°*  5e  +  °*  o625e3  -  0.0l30208e5  -  0. 77582465  X  10"2  e? 

-  6.  16930  X  10'3e9  -  4.  96735185  X  10'3en 

-  4.0929341  X 10‘3e13  -  3.44089913  X  10‘3e15  +  . 
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Table  C-l 


THE  F ^mp  AND  G|pq  FUNCTIONS  (Continued) 

G31q  =  1  +  2e2  +  3.  734375e4  +  5. 769097222e6  +  8.  096605089e8 
+  10. 6829796e 10  +  13. 50467544e12  +  16. 54355294e14 
+  19. 7851081e16  +  23 . 2 1740226e 1 8  +  26.  83038574e20  + 


G322  =  1*375e2  + 3.0625e4  + 5.099283854e6  +  7.427332899e8 

+  10. 01418813e10  +  12.  83625022e 12  +  1 5. 875406 17e14 

+  19. 1 17176 18e 18  +  22.  54963 923e 18  +  26. 16275785e20  + 

G30  2  =  °* 125eZ  +  °*  0208333e4  +  0. 0179036e6  +  0. 0127713e8 

+  0. 009776e10  +  0. 0077966e12  +  0. 0064085e14  + . 

G41  l  =  0.  5e  +  2. 0625e3  +  4. 794271e5  +  7.048177e7 

+  7. 4l6341e9  +  6.046672eJI  +  3.  999685e13 

+  2.208745e15  +  1.037865e1?  +  0.420094e19  + . 

G.ln  =  1  +  e2  +  4.  0625e4  +  6. 3125e6  +  6. 566406e8 
410 

+  5.257813e10  +  3.438477e12  +  1.  913086e14 

+  0. 923584e16  +  0.  389633e18  +  0. 1436l3e20  + . 

°421  =  2*  5e  +  8-4375e3  +  18.65885e5  +  26.  132813e7 

+  26.  069336e9  +  20.032104eU  +  12.  389781e13 

+  6.332082e15  +  2.7l6932e17  +  0. 986036e19  + . 

G40  2  *  0*5e2  -  0.3333333e4  -  0.125e6  +  0. 338542e8 
+  0.  3  14453e 10  +  0.09309<re12  -  0.048706e14 
-  0.062805e16  -  0.027145e18  -  . 
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Table  C-2 


VALUES  OF  LOW-ORDER  (FOURTH)  HARMONICS  OF  EARTH 
GRAVITATION  POTENTIAL  (WGS  72  MODEL) 


Zonal  Harmonics 
J2  =  1082.61579  E-6* 

J3  =  -2. 53881  E-6 
J4  =  -1.  65597  E-6 

Tessera l  Harmonics 

S22  =  _9'  0602  E"7 

531  =  2.8843  E-7 

532  =  02.2055  E-7 

533  =  1.  9611  E-7 
S42  =  1.4562  E-7 
S..  =  6. 7006  E-9 


C22  =  1.  5765  E-6 
C31  =  2.2181  E-6 
C32  =  3.  1196  E-7 
C33  =  9.8324  E-8 
C42  =  7.6894  E-8 
c44  =  *4. 0641  E-9 


